## Same Policy, Different Dreams: Explaining the Variation in Local Net-zero Policy Stringency in South Korea
## Authors: Yoonsoo Kim, Gyuhwan Kim, Inhwan Ko
## Date: Mar 18, 2025

#install.packages("pacman")
pacman::p_load(MASS, tidyverse, Gifi, readxl, sandwich, lmtest)

## 1. Create the dependent variable using Gifi (non-linear PCA)
data <- read.csv("korord_data.csv")

# 1) target year
data$x1_targetyear[data$targetYear==2050] <- 0
data$x1_targetyear[data$targetYear!=2050] <- 1
sum(is.na(data$x1_targetyear))

# 2) interim target
data$x2_interim[data$targetInterim=="no"] <- 0
data$x2_interim[data$targetInterim!="no"] <- 1
sum(is.na(data$x2_interim))

# 3) annual plan
data$x3_annualplan[data$targetPlanAnnual=="no"] <- 0
data$x3_annualplan[data$targetPlanAnnual=="yes"] <- 1
sum(is.na(data$x3_annualplan))

# 4) target monitor
data$x4_monitor[data$targetMonitor=="no"] <- 0
data$x4_monitor[data$targetMonitor=="annual"] <- 1
sum(is.na(data$x4_monitor))

# 5) committee - chair is mayor?
data$x5_chair[data$comChair=="no"] <- 0
data$x5_chair[data$comChair=="yes"] <- 1
sum(is.na(data$x5_chair))

# 6) committee - member is diverse?
data$x6_memdiv[data$comMemberDiv=="no"] <- 0
data$x6_memdiv[data$comMemberDiv=="yes"] <- 1
sum(is.na(data$x6_memdiv))

# 7) inventory
data$x7_inven[data$inventory=="no"] <- 0
data$x7_inven[data$inventory=="yes"] <- 1
sum(is.na(data$x7_inven))

# 8) carbon neutrality fund
data$x8_fund[data$fund=="no"] <- 0
data$x8_fund[data$fund=="yes"] <- 1
sum(is.na(data$x8_fund))

# 9) just transition
data$x9_justTransition[data$justTransition=="no"] <- 0 
data$x9_justTransition[data$justTransition=="yes"] <- 1
sum(is.na(data$x9_justTransition))

# 10) carbon neutrality center
data$x10_center[data$center=="no"] <- 0
data$x10_center[data$center=="yes"] <- 1
sum(is.na(data$x10_center))

# 11) carbon neutrality special administrator
data$x11_admin[data$admin=="no"] <- 0
data$x11_admin[data$admin=="yes"] <- 1
sum(is.na(data$x11_admin))

prindata <- data[,35:45] %>% as.data.frame()
fitord <- princals(prindata, ndim=3, levels="nominal")

summary(fitord, cutoff=.01)
ordscore <- fitord[["objectscores"]][,1]
data$ordscore <- ordscore

## variables
data$partisanship %>% unique()
data$ideo[data$partisanship=="minju" | data$partisanship=="jinbo"] <- "liberal"
data$ideo[data$partisanship=="none"] <- "neutral"
data$ideo[data$partisanship=="kukmin" | data$partisanship=="libertykor"] <- "conserv"
data$ideo <- as.factor(data$ideo)

data$renewpot <- 100*(data$renewgen/ (365*24)) / data$renewcap

## manufacture
manumine <- read_excel("manumine.xlsx")
data <- left_join(data, manumine, by=c("koraff","korname"))
data$manumine_grdp <- data$manumine / data$grdp2019 *100

## affiliated high-level governments' score

affs <- data$aff %>% unique()
affs_score <- data %>% 
  select(name, ordscore) %>% 
  filter(name %in% affs)
colnames(affs_score) <- c("aff","affscore")
affs_score <- affs_score[-18,]

data$affscore <- NA

for (i in 1:length(affs)) {
  output <- affs_score$affscore[affs_score$aff==affs[i]]
  input <- data$affscore[data$aff==affs[i]]
  outputs <- rep(output, length(input))
  data$affscore[data$aff==affs[i]] <- outputs
}

data$elecuse <- gsub(",","",data$elecuse)
data$elecuse <- as.numeric(data$elecuse)

data$grdp2019[data$korname=="세종시"] <- 11863000

## 2. Which factors are associated with ordscore? OLS
#rawdata <- data

## 2.1. both high-level and low-level
colnames(data)
# renewable energy share
model1 <- lm(ordscore ~ 
               manumine_grdp + renewshare +
                 relevel(ideo, ref="neutral") + 
                 tmn + 
                 log(density+1) + log(grdp2019+1) + autonomy,
               data)

summary(model1)
coeftest(model1, vcovHC=vcov(model1, type="HC2")) %>% round(3)
logLik(model1)

## 2.2. only low-level
data2 <- data %>% filter(type=="municipal")
# renewable energy share
model2 <- lm(ordscore ~ 
               manumine_grdp + renewshare +
                 relevel(ideo, ref="neutral") + 
                 tmn + 
                 log(density+1) + log(grdp2019+1) + autonomy +
                 affscore + as.factor(aff),
               data2)
summary(model2)
coeftest(model2, vcovHC=vcov(model2, type="HC2")) %>% round(3)
logLik(model2)

nrow(data2)

## 2.3. spatial model

library(spaMM)

data3 <- data2
colnames(data3)[23] <- "y"
colnames(data3)[24] <- "x"

model3 <- fitme(ordscore ~ manumine_grdp + renewshare + 
                  relevel(ideo, ref="neutral") + 
                  tmn + log(density+1) + log(grdp2019+1) + 
                  autonomy + affscore + as.factor(aff) +
                  Matern(1|x+y), 
                data3,
                family=gaussian(link = "identity"))
model3
model3_results <- summary(model3)
model3_results$beta_table %>% round(3)

### current status of net-zero policies 

tracker <- read_excel("tracker.xlsx") %>% 
  filter(actor_type=="City")

unique(tracker$end_target)
netzero_type <- c("Climate neutral", "Net zero", "Carbon neutral(ity)",
                  "GHG neutral(ity)", "Zero emissions", "1.5°C target",
                  "Zero carbon", "Net negative", "Climate positive")
tracker$netzero <- ifelse(tracker$end_target %in% netzero_type, 1, 0)
nrow(tracker[tracker$netzero==1,])/1185
nrow(tracker[tracker$netzero==0,])

nrow(tracker[tracker$end_target=="Absolute emissions target" |
               tracker$end_target=="Reduction v. BAU" |
               tracker$end_target=="Emissions intensity target",])

unique(tracker$end_target_status)

tracker %>% 
  filter(netzero==1) %>% 
  group_by(end_target_status) %>% 
  summarize(n=n())

unique(tracker$end_target_year) 

tracker %>% 
  filter(netzero==1) %>% 
  group_by(end_target_year) %>% 
  summarize(n=n())

tracker %>% 
  filter(netzero==1) %>% 
  group_by(interim_target_year) %>% 
  summarize(n=n())

tracker %>% 
  filter(netzero==1) %>%
  filter(interim_target_year!="2050" | !is.na(interim_target_year)) %>% 
  group_by(interim_target_target_emissions) %>% 
  summarize(n=n())

tracker %>% 
  filter(netzero==1) %>% 
  group_by(reporting_mechanism) %>% 
  summarize(n=n())

colnames(data)
data %>% group_by(admin) %>% summarize(n=n()) %>% 
  mutate(per = 100*n/191)

## Figure 2

top10 <- data %>% slice_max(ordscore, n=10)
bottom10 <- data %>% slice_min(ordscore, n=10)
topbot <- rbind(top10, bottom10)

topbot %>% 
  arrange(ordscore) %>% 
  ggplot(aes(x=reorder(name,ordscore), y=ordscore, fill=name)) +
  geom_bar(stat="identity") +
  geom_text(aes(label=round(ordscore, 3), hjust=0.5)) +
  coord_flip() +
  scale_fill_manual(
    values=c('gwangju'='#1B6DC5','nowongu'='#237EDF','dangjin'='#2392DF',
             'suwon'='#23AFDF', 'gyeonggi'='#23AFDF', 'gangwon'='#23AFDF',
             'seongnam'='#23AFDF', 'jeonbuk'='#23AFDF', 'seodaemungu'='#2BE3E6',
             'seoul'='#2BE3E6', 'jeonnam'='#2BE3E6', 'incheon'='#2BE3E6',
             'goyang'='#2BE3E6', 'busan'='#2BE3E6',
             
             'youseong'='#E6B62B','yongsangu'='#E6B62B','seongju'='#E6B62B',
             'namdonggu'='#E6B62B', 'jeongeup'='#E6B62B', 'goseong'='#E6B62B',
             'dobonggu'='#E6B62B', 'taeahn'='#E69A2B', 'muju'='#E68E2B',
             'gwacheon'='#E68E2B', 'yangpyeong'='#E6662B', 'bukgu'='#E6662B'
             )
  ) +
  theme_light(base_size=15) +
  ylab("Local net-zero policy index") +
  xlab("Localities") +
  theme(legend.position="none")

## descriptive statistics

round(mean(data$manumine_per),2); round(sd(data$manumine_per),2)
round(mean(data$renewshare),2); round(sd(data$renewshare),2)

unique(data$ideo)
data$lib <- ifelse(data$ideo=="liberal", 1, 0)
round(mean(data$lib),2); round(sd(data$lib),2)

round(mean(data$tmn),2); round(sd(data$tmn),2)
round(mean(data$density),2); round(sd(data$density),2)
round(mean(data$grdp2019),2); round(sd(data$grdp2019),2)
round(mean(data$autonomy),2); round(sd(data$autonomy),2)
round(mean(data$manumine_grdp),2); round(sd(data$manumine_grdp),2)
